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LONG-TERM  GOALS 

In  this  work,  we  will  employ  a  high-resolution  coupled  hydrodynamic-sediment  model  to  examine  the 
relative  importance  of  the  principal  mechanisms  controlling  the  morphodynamics  of  Skagit  Bay.  Using 
the  measurements  from  the  extensive  observation  program  supported  through  the  tidal  flats  DRI,  we 
will  examine  the  capability  of  a  state  of  the  art  coastal  ocean  model,  and  determine  what  future 
extensions  may  be  necessary  for  continued  discovery  in  this  field.  Through  extensive  grid  refinement 
efforts  and  available  high-fidelity  bathymetry,  a  better  understanding  of  the  mesh  resolution  required  to 
resolve  the  critical  processes  will  be  gained.  This  will  guide  future  application  of  this  class  of  model. 

OBJECTIVES 


In  this  project,  we  will  configure  an  advanced  coupled  hydrodynamic-sediment  model  for  simulation 
of  the  circulation  and  sediment  transport  in  Skagit  Bay.  The  model  will  resolve  the  range  of  required 
scales  from  the  open  boundary  in  Puget  Sound  (~  50  km)  to  the  channel  networks  on  the  flats  (-10-100 
m).  The  coupled  model  will  be  validated  using  available  measurements  to  detennine  the  capabilities 
and  needs  of  such  a  system  for  this  class  of  application.  Grid  convergence  studies  will  be  performed  to 
determine  the  necessary  mesh  resolution  required  to  resolve  the  dominant  processes.  We  will  employ 
the  calibrated  coupled  model  to  evaluate  the  relative  importance  and  influence  of  observed  external 
forcing  (fluvial,  tidal,  wind,  wind-wave  and  surface  heating)  on  sediment  dynamics  and  morphological 
change  of  the  inter-tidal  region  of  Skagit  Bay  over  a  range  of  time  scales  from  tidal  to  seasonal. 

APPROACH 

We  will  develop  and  apply  a  coupled  hydrodynamic-sediment  model  of  Skagit  Bay.  Due  to  the 
complexity  of  the  coastline  and  bathymetry  and  the  large  range  in  dynamical  scales  in  macrotidal 
estuaries,  the  unstructured-grid  coastal  ocean  model  FVCOM  was  selected.  FVCOM  is  a  Fortran90 
software  package  for  the  simulation  of  ocean  processes  in  coastal  regions  (Chen  et  ah,  2003,  2006). 

The  publicly  available  model  has  a  growing  user  base  and  has  been  used  for  a  wide  variety  of 
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applications,  including  work  in  Skagit  Bay  (Yang  and  Khangaonkar,  2008).  The  kernel  of  the  code 
computes  a  solution  of  the  hydrostatic  primitive  equations  on  unstructured  grids  using  a  finite-volume 
flux  formulation.  For  the  vertical  discretization,  a  generalized  terrain-following  coordinate  is 
employed.  The  model  is  fully  parallelized  using  a  Single  Program  Multiple  Data  (SPMD)  approach 
(Cowles,  2008).  FVCOM  will  be  coupled  with  the  Community  Sediment  Transport  Modeling  System 
(http://woodshole.er.usgs.gov/project-pages/sediment-transport/).  The  model  includes  transport  of 
both  the  suspended  load  and  bedload.  The  number  of  sediment  classes  is  flexible,  and  for  each  class, 
parameters  such  as  critical  shear  stress,  mean  diameter,  and  settling  velocity  must  be  defined.  Complex 
bed  dynamics  are  included  with  a  user-prescribed  number  of  layers  defined  by  the  layer  number, 
fractions  of  each  sediment  class,  an  age,  and  a  thickness.  Recently,  support  for  cohesive  sediment 
modeling  has  been  implemented  in  CSTMS  and  is  currently  undergoing  testing. 

Beginning  from  the  aforementioned  FVCOM  model  of  Skagit  Bay  developed  at  PNNL,  we  will  work 
to  modify  this  model  to  resolve  the  processes  that  are  the  focus  of  the  proposed  work.  This  work 
includes  the  initialization  and  parameterization,  calibration  of  the  sediment  model,  modification  of  the 
mesh  resolution  and  expansion  of  the  domain  outward  from  Skagit  Bay.  The  updated  model  will  be 
validated  using  measurements  obtained  through  observational  programs  in  the  tidal  flats  DRI. 

An  open  loop  mesh  adaptation  procedure  will  be  used  to  refine  the  initial  model  setup  to  make  sure  the 
mesh  has  the  required  grid  spacing  in  tidal  channels  and  flats  needed  to  resolve  the  processes  examined 
in  the  proposed  work.  This  will  ensure  that  we  have  a  grid  refined  solution  for  our  morphodynamic 
studies  to  constrain  the  influence  of  spatial  discretization  errors  in  the  model  results.  To  do  so,  we 
will  generate  a  series  of  coarser  meshes  from  the  initial  grid  using  the  method  of  Maximal  Independent 
Sets.  Solutions  for  the  baseline  case  will  be  obtained  on  the  coarser  meshes  and  Richardson 
extrapolation  of  the  bed  thickness  change  and  integrated  current  energy  across  the  mesh  levels  will  be 
utilized  to  determine  an  estimate  of  local  spatial  discretization  error.  The  mesh  will  be  adapted  locally 
by  inserting  and  deleting  points  based  on  this  heuristic.  This  process  will  be  repeated  until  suitable 
grid  convergence  is  achieved  for  the  given  forcing  and  available  bathymetry.  Using  this  method  we 
will  take  advantage  of  the  variable  resolution  allowed  by  the  unstructured  mesh,  and  rely  on  the 
dynamics  to  drive  the  background  length  scale. 

We  will  force  the  model  using  idealized  and  realistic  conditions  to  examine  the  contributions  of 
external  forcing  to  the  morphodynamics  of  the  tidal-flat.  In  these  process  studies  the  rates  of 
freshwater  discharge  and  fluvial  sediment  supply  will  be  varied,  along  with  the  tidal  phase,  and  the  sea 
level  to  examine  the  impact  on  net  sediment  transport  and  circulation  in  Skagit  Bay.  In  addition,  we 
will  examine  the  impact  of  resuspension  driven  by  episodic  wind  and  wave  events.  These  studies  will 
be  used  to  examine  how  the  intertidal  zone  adjusts  on  event  time  scales,  through  the  spring-neap  cycle, 
and  on  seasonal  time  scales. 

WORK  COMPLETED 

1.  Model  setup 

Following  acquisition  of  the  FVCOM  Skagit  Bay  model  from  PNNL,  the  model  has  been  been 
modified  to  make  it  suitable  for  answering  the  scientific  questions  of  the  DRI.  This  includes  an 
increase  in  the  resolution,  modification  of  the  domain,  and  improvement  of  the  bathymetry.  These  are 
outlined  in  more  detail  below. 
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Domain 


The  model  boundaries  have  been  extended  into  Juan  de  Fuca  strait  and  through  the  Saratoga  Passage  to 
the  south  end  of  Whidbey  Island  at  Sandy  Point  (Fig  1).  The  extensions  were  performed  to  ensure  that 
the  open  boundary  regions  outside  Skagit  Bay  contained  had  a  large  enough  volume  relative  to  the 
tidal  prism  to  reduce  the  influence  of  the  open  boundary  hydrography  and  sediment  concentration  on 
the  interior  model  solution.  The  model  is  forced  by  tides  at  three  open  boundaries:  Sandy  Pt,  Juan  de 
Fuca,  and  Swinomish. 


Figure  1:  Skagit  Domain  and  Bathymetry  (Google  Earth  overlay) 


Bathymetry 

The  bathymetry  consists  primarily  of  two  components,  the  Puget  Sound  Digital  Elevation  Map 
(PSDEM2005)  with  full  coverage  of  the  model  domain  and  the  SRSC  Fir  Island  Lidar  survey, 
primarily  constrained  to  the  upper  flats  along  Fir  Island.  This  data  was  projected  and  reinterpolated 
into  geographic  coordinates  by  Rich  Signed  (USGS  Woods  Hole)  and  is  available  for  download  and 
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interrogation  through  several  mechanisms  via  the  USGS  Thredds  catalog  (http://coast- 
enviro.er.usgs.gov/thredds/bathy_catalog.html,  Puget  Sound  DEM:  2005). 

On  the  South  Fork  flats,  the  bathymetry  is  still  influenced  by  artifacts  present  in  the  DEM.  Several 
attempts  were  made  to  eliminate  these  artificts  which  manifest  themselves  as  lines  of  discontinuity  in 
the  bathymetry  in  roughly  the  along-flat  direction  at  intervals  of  roughly  a  kilometer.  At  this  time 
these  have  principally  been  dealt  with  by  smoothing  in  the  cross-flat  direction  to  try  to  maintain  a 
reasonable  large-scale  cross-flat  slope  without  smoothing  channel  features.  Digital  elevation  in  the 
primary  river  beds  has  been  estimated  in  collaboration  with  Dave  Ralston  (WHOI).  Channel 
centerlines  for  the  North  Fork,  South  Fork,  and  Main  Branch  of  the  Skagit  River  up  to  Mt.  Vernon 
along  with  Steamboat  and  Freshwater  Slough  were  provided  to  D.  Ralston.  De-tided  depth  collected 
from  the  WHOI  survey  (Geyer  et  al.)  is  then  interpolated  on  the  tracks  and  applied  to  model  rivers 
under  the  assumption  of  zero  cross-river  gradients.  While  this  method  produced  reasonable  slopes, 
some  modification  at  the  North/South  fork  split  was  necessary  to  maintain  continuity.  Deepwater 
Slough  was  graded  in  the  same  manner,  although  at  this  time  there  is  limited  flux  through  the  Slough. 
Bathymetry  in  the  swinomish  is  modified  to  enforce  the  dredged  controlling  depth  of  6.8  feet. 

Forcing 

The  model  is  forced  using  freshwater  flux  from  the  Mt.  Vernon  stream  gauge  and  specified  sea  surface 
elevation  (at  present,  tidal  elevation)  at  the  open  boundaries.  Additional  experiments  will  use 
temporally-varying  wind,  wave-induced  bottom  stress,  and  sediment  load. 

Model  Versions 

Two  series  of  models  have  been  generated.  The  3.xx  series  models  are  generated  using  the  mesh 
generation  component  of  the  Surface  Modeling  System  (SMS)  software  which  uses  an  advancing  front 
technique  to  generate  Delaunay-conforming  grids.  Working  with  D.  Ralston  (WHOI)  these  models 
have  been  continually  refined  to  ensure  stability  in  the  integration  while  maintaining  the  level  of 
necessary  resolution.  Currently  two  versions  (3.15,  3.16)  are  being  used  by  D.  Ralston  for  hindcast 
simulations  of  Skagit  Bay.  A  summary  of  metrics  for  these  models  is  included  below.  The  4.xx  series 
of  models  are  generated  using  the  open  source  meshing  software  gmsh  (http://www.geuz.org/gmsh/  ). 
This  software  has  the  capability  of  performing  automated  adaptive  mesh  refinement.  The  mesh  can  be 
adapted  locally  by  specifying  a  background  length  scale  to  gmsh.  Version  4.3,  a  relatively  coarse  mesh 
(table  1)  has  been  used  for  calibrating  and  validating  the  tidal  dynamics. 

Table  1:  Specifications  of  Current  Skagit  Models 


Model 

Elements 

Layers 

Resolution 

[Flats/Mean] 

Time  Step  (s) 

iterations/ (c  ore -hour) 

skg3.15 

188K 

21 

10m  /  37m 

0.5 

150 

skg3.16 

1 12K 

21 

20m  /  50m 

1 

255 

skg4.3 

15K 

11 

100m  /  200m 

5 

7000  (baro tropic) 
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2.  Non-Cohesive  dynamics  from  CSTMS  implemented  into  FVCOM 

Dynamics  for  cohesive  sediment  transport  have  been  added  to  the  sediment  model  in  FVCOM.  These 
have  been  extracted  primarily  from  the  implementation  of  the  Community  Sediment  Transport 
Modelling  System  (CSTMS)  currently  available  in  the  CSTMS  branch  of  the  Regional  Ocean  Model 
System  (ROMS).  The  principal  modification  to  the  existing  FVCOM  sediment  model  was  the 
treatment  of  the  bed.  For  cohesive  substrates,  the  critical  shear  stress  is  a  property  of  the  bed  which 
evolves  due  to  consolidation,  swelling,  and  bioturbation. 

To  test  the  updated  FVCOM  sediment  model  we  have  developed  a  wrapper  to  the  module  to  drive  1-D 
solutions  of  sediment  in  the  water  column.  Time-dependent  physical  fields  (eddy  diffusivity,  bottom 
stress,  and  depth)  are  generated  by  the  General  Ocean  Turbulence  Model  (GOTM)  and  are  used  to 
drive  the  model.  This  enables  controlled  testing  of  dynamics  and  sediment  algorithms  for  problems 
that  vary  slowly  in  the  horizontal.  With  C.  Sherwood  (USGS,  Woods  Hole)  we  are  in  the  process  of 
developing  a  standardized  set  of  1-D  test  cases  for  the  purpose  of  debugging  as  well  as  testing  the 
CSTMS  model  dynamics  in  FVCOM.  This  cases  will  explore  steady  state  forcing  as  well  as  tidal 
(mixed  tides)  and  event-scale.  With  guidance  from  C.  Sherwood  we  are  also  examining  several 
possibilities  for  inclusion  of  flocculation  dynamics  within  the  CSTMS  framework,  focusing  on 
methods  that  account  for  redistribution  of  sediment  across  classes  associated  with  floe  breakup  and 
formation  (Xu  et  ah,  2008). 

RESULTS 

1.  Tidal  Modeling 

Considering  the  depth  of  water  and  the  relative  protection  and  short  fetch  of  the  Bay,  the  tides 
represent  the  dominant  external  forcing  over  the  broad  scale  of  the  flats  and  thus  must  be  simulated  as 
accurately  as  possible.  However  due  to  the  extensive  flooding  and  drying  of  the  flats,  the  lateral 
mixing  associated  with  the  extreme  currents  of  Deception  Pass  and  the  diurnal  inequality,  achieving 
high  level  of  tidal  model  skill  is  challenging.  A  coarse  model  has  been  used  for  testing  the  model’s 
ability  to  reproduce  the  tidal  harmonics  at  stations  around  the  Bay.  This  model  runs  quickly  (Table  1) 
and  thus  is  useful  for  tidal  calibration.  This  model  can  be  used  to  inform  the  forcing  and  setup  of  the 
high-resolution  models  which  are  too  large  to  be  practical  for  anything  but  production  runs.  Two  key 
with  the  tidal  forcing  are  the  selection  of  open  boundary  forcing  at  the  Juan  de  Fuca  open  boundary 
where  there  are  no  available  harmonics  and  the  treatment  of  Deception  Pass.  Tidal  harmonics  are 
compared  at  9  stations  in  the  domain:  Deception  Pass  Park,  Holly  Farms  Harbor,  GreenBank, 
Coupeville,  Crescent  Harbor,  Sneeoosh  Point,  Ala  Spit,  Corney  Bay,  and  Yokeko  Point.  For  the  M2 
amplitude  and  phase  (dominant  constituent)  the  average  amplitude  error  is  3.3  cm  and  the  average 
phase  error  is  4.5°.  The  primary  source  of  model-observation  error  is  the  upper  Skagit  where  the 
model  leads  the  observations  in  the  M2Constituent  by  approximately  6°.  We  are  continuing  work  on 
parameterization  of  the  bottom  friction  in  Deception  Pass  to  improve  the  tidal  response  in  the  upper 
Skagit.  Spatially  distribution  of  model-computed  M2  amplitude  (m)  and  phase  (°G)  are  shown  in 
Figure  2. 
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Figure  2:  Model-Computed  M2  Amplitude  (m)  [left]  and  phase  ( °G)  [right]  in 

water  greater  than  2-m  depth 


2.  FVCOM-CSTMS  testing 

With  guidance  from  C.  Sherwood  (USGS),  a  suite  of  1-D  test  cases  for  the  FVCOM-CSTMS  model 
are  being  developed.  These  cases  will  be  used  to  debug  and  to  test  the  sediment  dynamics  module  of 
FVCOM  using  a  range  of  forcing  scenarios.  The  test  cases  allow  evaluation  of  parameter  sensitivity 
and  comparison  with  RO. 


Table  3:  FVCOM-CSTMS  1-D  Test  Cases 

Test  Case  Open  Channel  Tidal  Forcing  Event  Forcing  Settling  Chamber 


Source 

Warner  et  al, 

2008 

Sherwood  et  al, 

2008  ppt 

Sherwood  et  al, 

2008  ppt. 

Sherwood  et  al., 

2008  ppt. 

Tests 

Steady-state 

equilibrium, 

turbulence 

models 

Time-varying 
erosion/deposition, 
bed  composition 

Dynamics  across 
two  events 
including  bed 
dynamics 

Quiescent  settling, 
steady  state  bed 
composition 

Physical 

Forcing 

Constant 

Pressure 

Gradient 

Harmonic  Pressure 
Gradient 

Time-varying 
Surface  Stress 

None 

Sed  Classes 

Sand 

Mud/Mud 

Mud/Mud 

Mud/Mud 

6 


u  (m/s) 


Figure  3:  Comparison  of  FVCOM-CSTMS  (red)  and  ROMS-CSTMS  (black)  for 

the  open  channel  sediment  test  case 


IMPACT/APPLICATIONS 

A  key  component  of  this  work  is  the  development  of  an  unstructured  grid  coupled  hydrodynamic- 
sediment  model  for  the  study  of  tidal  flats.  The  validation  efforts  will  draw  on  the  intensive 
observation  program  supported  by  this  DRI  and  will  help  to  make  conclusions  about  the  potential  of 
such  a  model  as  well  as  define  future  research  needs  in  terms  of  development  or  need  for  additional 
data.  In  concert  with  the  scientific  findings  of  other  teams,  the  coupled  system  will  provide  a  platform 
by  which  the  influence  of  external  forcing  components  on  tidal  flat  sediment  processes  can  be  isolated 
and  elucidated. 

RELATED  PROJECTS 

In  this  work  we  work  closely  with  other  investigators  participating  in  the  ONR  tidal  flats  DRI.  The 
key  collaborators  are  C.  Sherwood  and  R.  Signed  (USGS,  Woods  Hole)  who  will  be  assisting  with  the 
development,  implementation,  and  validation  of  CSTMS  within  FVCOM  as  well  as  processing  of 
bathymetry  for  the  model  domain.  We  are  also  working  closely  with  D.  Ralston  (WHOI)  and  J. 
Lerczak  (OSU)  in  the  model  development  and  validation. 
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